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Abstract 

Grand canonical Monte Carlo simulations are used to explore the metastable fluid-fluid coexistence 
curve of the modified Lennard-Jones model of globular proteins of ten Wolde and Frenkel (Science, 
277, 1975 (1997)). Using both mixed-field finite-size scaling and histogram reweighting methods, 
the joint distribution of density and energy fluctuations is analyzed at coexistence to accurately 
determine the critical-point parameters. The subcritical coexistence region is explored using the 
recently developed hyper-parallel tempering Monte Carlo simulation method along with histogram 
reweighting to obtain the density distributions. The phase diagram for the metastable fluid- fluid 
coexistence curve is calculated in close proximity to the critical point, a region previously unattained 
by simulation. 
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I. INTRODUCTION 



Many pathological diseases, including sickle cell anemiaii and genetic cataracts^, are 
known to be caused by the crystallization of certain globular proteins. The importance of 
proteins is further exemplified by recent advances in genome sequencing, revealing that as 
much as ninety-eight percent of DNA may be involved in the regulation of their production. 
Exploring protein structure and activities (proteomics) is a growing research field and should 
help in our understanding of health and disease on a molecular basis^. Advances in decod- 
ing genomes, however, have far and away surpassed those in the determination of protein 
structure. The growth of high quality protein crystals from solution is important in deter- 
mining structure and is known to depend sensitively on the initial conditions of the solution. 
Unfortunately, knowledge of the initial conditions necessary for optimal crystallization has 
not come easily^. 

Significant progress in understanding the relationship of the initial conditions to the 
crystal nucleation rates for globular proteins has been made in recent years. It was demon- 
strated in the theoretical work by Gast, Hall, and Russell that the characteristics of the 
phase diagram of colloids depend sensitively on the range of attraction between the col- 
loidal particles. It was found that for a colloid-colloid attractive interaction of very short 
range (less than some thirty percent of the colloidal diameter), there are stable fluid and 
crystal phases and a metastable fiuid-fiuid phase located below the liquidus line. Other 
studies^ have found similar results. This has also been demonstrated in both experimenl!^ 
and simulation^. Rosenbaum, Zamora, and Zukoski^ linked the experimental observations of 
George and Wilsoi>i^ with those of colloids. They found that the narrow range in the second 
virial coefficient for which globular proteins crystallize map onto an effective temperature 
range for colloidal systems. The phase diagrams of the two systems were analogous, ten 
Wolde and FrenkeUi then calculated the phase diagram and nucleation rate for a modified 
Lennard- Jones model of globular proteins, whose range of attractive interaction was small 
compared with the protein diameter. In this seminal work they showed that the nucleation 
rate increased by many orders of magnitude in the vicinity of the critical point, suggesting a 
direct route to effective crystallization. Therefore, accurate knowledge of the region around 
the critical point provides important information regarding crystallization. 

ten Wolde and Frankel— studied a modified Lennard- Jones (ML J) pairwise interaction 
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V(r) is shown in Fig. 1, where a denotes the hard-core diameter of the particles, r is the 
interparticle distance, and e is the well depth. The width of the attractive well can be ad- 
justed by varying the parameter a; it was tuned in such a way that the so-called stickiness 
parameter was equal to that produced by the hard-core Yukawai^ potential for k = 7a~^ 
at the metastable liquid-vapor critical point, where is a measure of the range of the 
attractive part of the potential. In the Yukawa model, the phase diagram for k = 7a^^ was 
found to be equivalent to that of globular proteins and maps onto those determined experi- 
mentally^. An advantage of the MLJ potential is that it lends itself naturally to both Monte 
Carlo and molecular dynamic simulations. Of particular interest both theoretically^»ii*i^ and 
experiment allyi^*^^*^ is the metastable fiuid-fiuid curve of the phase diagram, for reasons 
noted above, ten Wolde and Frenkel^i determined the phase diagram for the above model 
using the Gibbs ensemble Monte Carlo (GEMC) methoc^^ where two coexisting phases are 
separated into two physically detached but thermodynamically connected boxes, the vol- 
umes of which are allowed to fluctuate under constant pressure. In the neighborhood of the 
critical point, however, GEMC cannot be relied upon to provide accurate estimates of the 
coexistence curve parametersii^iSiiS. This is evident by considering the metastable region of 
the phase diagranJ^. The error bars in the data points grow larger as the critical point is 
approached; there are no data points in the immediate vicinity of the critical point. The 
purpose of this paper, then, is to accurately determine the critical point of the phase diagram 
of ten Wolde and Frenkel using finite size scaling techniques adapted for simple fluids by 
Bruce and Wilding and to accurately determine the corresponding subcritical region using 
the hyper-parallel tempering method. 

II. MODEL 
A. Theory 

To study the critical region of the phase diagram, we use the Bruce and Wilding^S 
finite-size scaling method, along with histogram-reweighting^i techniques, in conjunction 
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FIG. 1: The potential in Eq. (1) is shown as the sohd hne. For comparison, the traditional 
Lennard-Jones potential is shown as the dashed line. The width of the MLJ was chosen so as to 
agree with the phase diagrams of experimentally determined globular proteins. 

with grand canonical Monte Carlo (GCMC) simulations^^. We write the reduced chemical 
potential as fi = fi/^ksT) and the reduced well depth as u = e/^ksT), with T the temper- 
ature of the system. In what follows, we denote fi as fi for simplicity. The critical point is 
identified by the critical values of the reduced chemical potential /ic and well depth Uc- The 
relevant scaling fields comprise linear combinations of /i and uj: 

T = Uc- UJ + S{fl - flc) (2) 




MLJ 
U 



h = n- Hc + riuJc-oj), (3) 

with r and h the thermal and ordering scaling fields, respectively. The parameters r and s 
are system-specific parameters, controlling the degree of mixing^ and vanish identically when 
the Ising symmetry is present. Conjugate to these two fields are the scaling operators 

M = j^^[p-su] (4) 

with M. and £ the ordering and energy-like operators, respectively. The (dimensionless) 
number density and energy density are defined by p = L~'^Na'^ and hy u = U / iVe), respec- 
tively, with U the total energy of the system, V = L'^ the volume of the system, and e the 
well depth of the potential energy. 
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We make the usual finite-size scaling ansatz^: 



Pl{M,£) ^ PmA^m ^^^4 (6) 



where 



and 



= aeL'/", Am = aML^'^", ^m^^m = ^^4 = L\ (7) 



5M = M- {M)c S£ = £- {£)„ (8) 



where the parameters and as are non-universal scaling factors, (3 and u are the criti- 
cal exponents for the coexistence curve and the correlation lengtb^^, respectively, and the 
subscript c in the above equations denotes that the averages are to be taken at criticality. 
From our simulations using GCMC, we obtain the joint probability distribution of p and u 
at a particular point in parameter space of inverse temperature ( T) (3 and reduced chemical 
potential /i from which we obtain the joint probability distribution of mixed operators via 

PL{M,£) = {l-sr)PL{p,u). (9) 
Integrating out the energy-like dependence from the latter distribution gives 

Pl{M) = j d£PL{M,£). (10) 

Using the fact ^^i^^i^^i^'^i^^i^^ that the critical behavior of this distribution is in the same 
universality class as the Ising model, we match the above probability distribution to that of 
the universal fixed point function 

P*m{x) = j PX,Ax,y)dy, (11) 

which is known from an independent study^S. The critical point of the fluid can thus be 
estimated by tuning the temperature T, chemical potential p, and field-mixing parameter s 
such that Pl(7V1) collapses onto -Pj^(x). 

To save computer-time and to facilitate the matching process, we employ histogram- 
reweighting^i in lieu of performing tedious simulations. Once the joint probability distri- 
bution at a particular point in parameter space (/?, p) is obtained in one simulation run, 
information around its neighboring points in parameter space {(3', p') is extracted using 



exp[(/i- - fi)pV - {(3' - /3)^y]Pf '")(p, u) 
E exp[(/i' - ^,)pV - (/?' - fi)uV\Pt'^ (p, u) ■ 
Histogram-reweighting provides an accurate determination of the new probability distribu- 
tion and is well suited for Monte Carlo studies. Its limitations are discussed elsewhere^. 

In the subcritical region, we employ hyper-parallel tempering^ Monte Carlo (HPTMC) to 
sample the joint probability distribution. Below the critical point, density fluctuations are no 
longer large enough for the joint probability distribution to be accessible by standard Monte 
Carlo simulations. Thus, a free energy barrier exists below Tc which needs to be overcome 
if one is to sample both coexisting phases. HPTMC allows one to effectively tunnel through 
this barrier by swapping particle configurations between different simulations (replicas) at 
different state points. Other techniques^ sample this region by artificially lowering this 
barrier. 

In the grand canonical ensemble, the partition function can be written as 

Z = ^ ^[x)exv{-fiU{x) + [iN{x% (13) 

where x denotes the state of the system, fi(x) is the density of states, U{x) is the total 
energy of state x, and A^(x) is the number of particles in state x, and all other variables are 
defined as above. In accord with practice^, we consider a composite ensemble consisting of 
M non-interacting replicas, each at a different set of state points. The partition function of 
the composite ensemble is specified by 

M 

Zc = l[Zi, (14) 

i=l 

where the complete state of the composite ensemble is given by 

X = {xi,X2, ...,xm), (15) 

with Xi denoting the state of the i^^ replica. The unnormalized probability density of the 
composite state x is given by 

M 

p(x) = J] exp[-pU{xi) + fiN{xi)]. (16) 

i=l 
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To sample configurations from tlie composite ensemble, a Markov chain is constructed to 
generate configurations according to the limiting function in Eq. (16). In the Markov chain, 
two types of trial moves are employed: 1) within each replica, insertion/deletion trial moves 
are attempted according to standard Monte Carlo as adapted for use in the grand canonical 
ensembl©2fi, and 2) Configuration swaps are attempted between pairs of replicas i and i+1 
such that 



new _ old new _ old f^7^ 

To enforce a detailed-balance condition, the pair of replicas that are attempted to be 
swapped are chosen at random, and the trial swap is accepted with the probability 



Paccixi ^ Xi+i) = min[l, exp{Af3AU - AfiAN)], (18) 

where A/3 = Pi+i-^i, AU = U{xi+i) -U{xi), A/i = /ij+i -/ij, and AA^ = N{xi+i) - N{xi). 
Once joint probability distributions are obtained in this way, histogram reweighting is applied 
to obtain coexistence according to Eq. (12) and an " equal- weight"«2i construction: 



P^p^^'Ip) = 0.5. (19) 

It should be recognized that the average density (p) in the upper limit of the integral in Eq. 
(19) is itself a function of temperature and chemical potential and can therefore be obtained 
from the first moment of the reweighted histogram. 



B. Computational Details 

We studied a system of N particles contained in a three-dimensional, periodic cubic 
simulation cell having a volume V = L'^. Two particles separated by a distance r interact 
via the modified Lennard- Jones (MLJ) pair potential given in Eq. (1), where e and a 
denote the energy and length scales, respectively. The total energy, f/, is obtained by 
summing over all distinct pairs of particles. We employ a truncated, unshifted version of 
Eq. (1) using a cutoff radius Tc = 2.0 a, in accorcP^ with ten Wolde and Frenkel. The 
simulations were performed on system sizes of L = 6 a, 7 cr, 8 cr, and 10 a, implemented using 
the Metropolis Monte Carlo method as adapted for use in the grand canonical ensemble^S 
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with a constant volume L"^, inverse temperature (T) /3, and reduced chemical potential /i. 
The thermodynamic potential needed is actually fi* = fi — 3 ln[A/a] (where A is the thermal 
deBroglie wavelength)^. As in other implementations^, only particle insertion and deletion 
steps were employed, particle displacements being realized within the cell from a succession 
of such steps. In our simulations, equilibrium times used were approximately two million 
steps and production runs ranged from 500 million steps for the smaller system sizes to one 
billion steps for the higher system sizes. Such a high number of steps was needed in order 
to obtain smooth distributions at the low temperatures we are studying. We emphasize 
that the true equilibrium state of the system is that of fluid and solid coexistence. GCMC, 
however, is limited to 'low' densities and, thus, the solid region of the phase diagram is 
unaccessible. 

Using a previous estimate of the critical temperature^^ for this model, we first attempted 
to locate the critical point by tuning the reduced chemical potential fi* until the density 
distribution exhibited a double peaked structure. Once obtained, longer runs were per- 
formed to accumulate better statistics. To obtain two-phase coexistence and appropriate 
values of the field-mixing parameter s, we adopted the criterion^S that the order distribution 
Pl(A^) = / d£ Pl{M.,£) be symmetric m M. — {M). This criterion is the counterpart of 
the coexistence symmetry condition for the Ising model magnetization distribution. Having 
obtained in this manner a two-phase distribution near the critical point, we then matched the 
order-operator distribution Pl{M.) to that of the universal fixed point distribution P^[x). 
Employing histogram-reweighting, we tuned the chemical potential /i*, temperature T, and 
mixed-field parameter s until our distribution Pl{M.) collapsed onto that of the fixed point 
distribution. Once that was attained, we then attempted to match Pl{£) to the correspond- 
ing energy-like fixed point distribution P^iy) by tuning the field-mixing parameter r. This 
procedure was repeated for the various system sizes studied. 

In the subcritical region, we used six replicas for the L = 7 a system to obtain joint 
probability distributions. Small runs of approximately five million steps were first performed 
to optimize the choice of /i* for each replica. At low temperatures, fi* was chosen so that the 
density distribution was biased toward a high-density peak; conversely, at high temperatures, 
fi* was chosen so that low-density peaks were favored. Such a choice of fi* allowed high- 
density configurations to be "melted" when passed to high temperatures, allowing for a 
fuller exploration of configuration spaced. Also, we chose these values so that a high swap 
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frequency between replicas was realized. On average, replicas were swapped twenty-five 
percent of the time. Once these criteria were satisfied, Monte Carlo runs were extended to 
approximately 150 million steps to acquire good statistics. Histogram-reweighting was then 
applied to the resulting joint probability distributions to obtain coexistence in the subcritical 
region in accordance with Eq. (12) and Eq. (19). 

III. RESULTS 

The resulting density distributions obtained in the prescribed manner at the size- 
dependent critical point Tc(L) and yUc(-^) shown in Fig. 2. As is evident, the distributions 
become narrower with increasing system size L, approaching the limiting form of the fixed 
point distribution P^{x). The corresponding energy- density distributions are shown in Fig. 
3. 

Both distributions show an asymmetry, due to field- mixing effects, with that of Pl{u) 
being much more pronounced. These effects die off with increasing L, so that the limiting 
forms of both Pl(p) and Pl{u) match the fixed point distribution P^(x). As noted, this 
limiting form is easily recognizable for the density distributions. However, the limiting form 
of the energy density distributions does not follow this pattern. This difference is attributable 
to the coupling that occurs for asymmetric systems between the ordering operator and 
energy-like operator fluctuations. Those of A4 dominate at large L, while those in S do 
no t^^i^^ . Thus, a 'background' effect perturbs the energy- density distributions. 

Our estimates for the fixed point distribution are shown in Fig. 4. As is seen from the 
figure, the agreement is good for all system sizes studied. This matching alone allows us to 
accurately determine the critical point Tc(L) and yUc(-^) ^o also obtain good estimates 
of the field-mixing parameter s. Though matching of Pl{^) to P^iy) should also give good 
estimates of the critical point, and the field-mixing parameter r, fluctuations in the energy- 
like operator S are relatively weak and therefore do not allow for good matching^. It is 
unknown presently how to remove this background effect for fluid-systems, though it has 
been removed for the Ising model^. Nevertheless, we can still obtain a rough estimate for 
the field-mixing parameter r by observing that the distributions Pl{£) should have a shape 
similar to P^{y). For completeness, we include these curves along with the energy- fixed 
point function in Fig. 5. 
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FIG. 2: Density distributions at Tc{L) and f^*{L) for the system sizes L = 6a, 8a, and Wa. For 
clarity, the distribution corresponding to the L = 7 system size is not shown. 
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FIG. 3: Energy distributions at Tc{L) and /U*(L) for the system sizes L = 6a, 7a, 8a, and 10a. 

From estimates of the finite critical-point temperatures Tc{L), we can estimate the critical- 
point for the infinite- volume system. Since contributions to Pl{M.) from finite values of r 
grow with system size like tL^^'^, the matching condition results in a deviation from the true 
critical point Tc(cxd) as 

T,(oo)-T,(L)ocL-(^+l)/^ (20) 

where 6 is the universal correction to scaling exponent^^ and u is the critical exponent for 
the correlation length^. In Fig. 6, we plot the apparent critical temperature Tc{L) as a 
function of L~^^~^^^^'^ . By extrapolating to infinite volume, we arrive at an estimate for the 
true critical temperature. Similar arguments in accounting for field-mixing effects^ apply 
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FIG. 4: The measured form of the ordering operator distribution Pl{M) for the system sizes 
studied. Shown for comparison is the universal fixed-point ordering operator distribution Pj^(x) 
(sohd hne). The data have been expressed in terms of the scaled variable 
where aj^ was chosen such that the distributions have unit variance. 
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FIG. 5: The measured form of the energy-like operator distribution Pl{£) for the system sizes 
studied. Shown for comparison is the universal fixed-point ordering operator distribution P^iy) 
(solid line). The data have been expressed in terms of the scaled variable y = a^^L'^~^f'^{£ — £c), 
where a^^ was chosen such that the distributions have unit variance. The lines are a guide to the 
eye. 

in obtaining the true critical density via 

(p),(L)-(p),(oo)(xL-('^-V'^), (21) 
where d is the dimensionality of the system. In Fig. 7, {p)c{L) is plotted as a function of 
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FIG. 6: The apparent reduced critical temperature plotted function of L-^^+i)/*^, with 9 = 
0.54— and u = 0.629^-. The extrapolation to infinite volume yields the estimate Tc = .4145(5). 



TABLE I: Summary of results for location of the critical point 

L UL) {p),iL) s i__ 

6 .4120(3) -2.939(1) .312(2) -0.120(3) -0.65(5) 

7 .4130(2) -2.932(2) .312(2) -0.115(2) -0.65(5) 

8 .4135(3) -2.930(2) .318(4) -0.110(1) -0.65(5) 
10 .4136(6) -2.932(3) .315(5) -0.110(1) -0.65(5) 



L ^1^^ to obtain an estimate of the true critical density. We summarize our findings in 
Table I. 

Using HPTMC, we were able to obtain joint distribution functions corresponding to 
points in the subcritical region of the phase diagram of ten Wolde and Frenkel. The choices 
we employed for the values of the temperature T and chemical potential /i* of each replica 
are displayed in Table II, along with the coexisting chemical potentials and coexisting 
densities and pi, obtained after reweighting. The corresponding density distributions 
in this region are shown in Fig. 8, where the temperatures are expressed in terms of the 
infinite-volume critical temperature. It can be seen that for temperatures even as high as 
T = .965 Tc, the system has a vanishingly small probability of visiting a state between the 
two coexisting phases. This implies, as stated previously, that one cannot expect standard 
Monte Carlo simulations to obtain these subcritical distributions. 

In Fig. 9, we plot the temperatures used in Table II as a function of density, where 
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FIG. 7: The measured average density {p)c{L) at the L-dependent critical point expressed as a 
function of L~'^'^~^^^'^ . Extrapolation to infinite volume yields the estimate pc = .3206(4). 
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FIG. 8: The coexistence density distributions for the L = 7 a system size for a range of subcritical 
temperatures, obtained using HPTMC as described in the text. 

we have calculated the average densities for each of the two phases from the distributions 
shown in Fig. 8. These new results are consistent with the earlier results for the metastable 
coexistence curve of the phase diagran>ii and extend much closer to the critical point. 
We also show our attempt to fit this subcritical density data to a power-law of the form 
p± — Pc = A\T — Tc\ ± B \T — Tc\^ , where and pc are are our extrapolated values as L — oo 
and f3 = .3258 is the Ising exponent. Although the data are in reasonable agreement with 
this fit, note that we have not accounted for possible finite-size effects, preventing this from 
being a definitive test. Fig. 10 shows the corresponding chemical potential as a function of 
temperature, which obeys a linear relationship in this region. Experimental^ investigations 



13 



TABLE II: Values of T and /i* used in HPTMC. Also shown are the reweighted chemical potentials 
used in obtaining the coexisting densities and pi. 



Replica T p* p*j^ py pi 

1 .400 -3.110 -3.137 .120 .531 

2 .404 -3.040 -3.071 .140 .510 

3 .406 -3.010 -3.041 .149 .495 

4 .409 -2.970 -2.995 .170 .473 

5 .4115 -2.947 -2.957 .186 .444 

6 .4123 -2.941 -2.946 .192 .441 



0.44- 



0.42- 




FIG. 9: Temperature plotted vs. density, as obtained from HPTMC simulations for the L = 7a 
system size shown with the critical point rc(oo) and the corresponding critical density pc as obtained 
by Eqs. (20) and (21). Also shown (sohd line) is a fit to the data, with /3 = .3258 i^. 

on the bovine lens protein, 7// — crystallin, have been performed near the critical point. 
They report values for the critical isothermal compressibility 7 = 1.21 ±0.05 and the critical 
correlation length i> = 0.68 ± 0.1. Both results are compatible to three-dimensional Ising 
model values. 
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FIG. 10: The chemical potential vs. T for the system size L = 7 a for the coexisting subcritical 
region obtained using HPTMC as described in the text. Errors do not exceed symbol size. 

IV. CONCLUSION 

In this work, we have employed mixed-field finite-size scaling techniques and histogram 
extrapolation methods to obtain accurate estimates for the critical-point parameters of the 
truncated, unshifted MLJ fluid with Tc = 2.0 cr. Our measurements allow us to pinpoint 
the critical-point parameters to within high accuracy. A previous estimate put the critical 
point at Tc = .42021, slightly higher than our estimate of Tc{oo) = .4145(5). In the near 
subcritical region, we have explored the phase diagram and obtained data in a region where 
no prior estimates were available. HPTMC proved to be an efficient means toward this end. 
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